%gridding with morpholigical method, preference:automatic analysis of DNA microarray images using mathematical
%p:input image , e.g. GSM16101_CH1-1.tif
function para=morphspotgrid1(p)
para=[];
st=cputime;
A=imread(p);
vv=0.5;
A=double(A);
A=A./256;
[h,w,d]=size(A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%horizontal projection
for i=1:h
    m=0;
    for j=1:w
          m=m+A(i,j);
    end
    Hy(i)=m/w;
end
%mean
hym=sum(Hy)/h;
%reduce mean
Hy1=Hy-hym;
%morphological reconstruction
Hy2=imreconstruct(Hy1,Hy);
clear Hy1;
%residue signal
Hy3=Hy-Hy2;
clear Hy2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%thresholding
up=vv*sum(Hy3)/h;
clear Hy;
%transfer the projection signal to binary signal
for k=1:h
    if Hy3(k)>up
        Hy4(k)=0;
    else
        Hy4(k)=255;
    end
end
clear Hy3;
%figure,plot(1:h,Hy4);
%vertical projection
for i=1:w
    n=0;
    for j=1:h
        n=n+A(j,i);
    end
    Vx(i)=n/h;
end
%mean
vxm=sum(Vx)/w;
%reduce mean
Vx1=Vx-vxm;
%morphological reconstruction
Vx2=imreconstruct(Vx1,Vx);
clear Vx1;
%redidue signal
Vx3=Vx-Vx2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%thresholding
up1=0.5*sum(Vx3)/w;
clear Vx2;
clear Vx;
%transfer the projection signal to binary signal
for k=1:w
    if Vx3(k)>up1
        Vx4(k)=0;
    else
        Vx4(k)=255;
    end
end
clear Vx3;
%figure,plot(1:w,Vx4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find coordinate y
hf=0;%flag
flag=0;
for i=1:h-2 %find the original pixel
    if flag==0
        if Hy4(i)==255 && Hy4(i+1)==255 && Hy4(i+2)==255
         hf=hf+1;
         flag=1;
         hspotcoff(hf)=i;
         t=0;
       end
    end
    if flag==1%find the terminal pixel
        t=t+1;
       if Hy4(i+1)==0 &&t>4
           hf=hf+1;
           hspotcoff(hf)=i;
           flag=0;
       end
    end
end
disp(['the original coordinate y��',num2str(hf)]);
clear Hy4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%compute the middle of two grid lines
disp('compute the middle to get final horizontal grid lines');
hspotcoff1(1)=hspotcoff(2);
for i=3:2:hf-1
    hspotcoff1((i+1)/2)=round((hspotcoff(i)+hspotcoff(i+1))/2);
end
hspotcoff1(round((hf-1)/2+1))=hspotcoff(hf);
clear hspotcoff;
%figure ,plot(1:w,Vx1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find coordinate x
vf=0;%flag
flag=0;
for i=1:w-2
    if flag==0
        if Vx4(i)==255 &&Vx4(i+1)==255 &&Vx4(i+2)==255
         vf=vf+1;
         flag=1;
        vspotcoff(vf)=i;
        t=0;
       end
    end
    if flag==1
        t=t+1;
       if Vx4(i+1)==0 &&t>4
           vf=vf+1;
           vspotcoff(vf)=i;
           flag=0;
       end
    end
end
disp(['original coordinate x��',num2str(vf)]);
clear Vx4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ompute the middle of two grid lines
disp('compute the middle to get final horizontal grid lines');
vspotcoff1(1)=vspotcoff(2);
for i=3:2:vf-1
    vspotcoff1((i+1)/2)=round((vspotcoff(i)+vspotcoff(i+1))/2);
end
vspotcoff1(round((vf-1)/2+1))=vspotcoff(vf);
hf=hf-1;
vf=vf-1;
clear vspotcoff;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
et=cputime-st;
disp(['gridding time:',num2str(et)]);
%draw the original grid line
[s,hf1]=size(hspotcoff1);
[s,vf1]=size(vspotcoff1);
up
up1
mm=5;
hspotcoff1(1)=hspotcoff1(1)-mm;
if hspotcoff1(1)<0
    hspotcoff1(1)=2;
end
hspotcoff1(hf1)=hspotcoff1(hf1)+mm;
if hspotcoff1(hf1)>h
    hspotcoff1(hf1)=h;
end
vspotcoff1(1)=vspotcoff1(1)-mm;
if vspotcoff1(1)<0
   vspotcoff1(1)=2;
end
vspotcoff1(vf1)=vspotcoff1(vf1)+mm;
if vspotcoff1(vf1)>w
    vspotcoff1(vf1)=w;
end
disp(['drar grid line','coordinate y��',num2str(hf1),',coordinate x��',num2str(vf1)]);
%Drawgrid(p,h,w,hspotcoff1,vspotcoff1,hf1,vf1);
%judge the gridding if correct or not, 1 is correct,0 is wrong
correct=iscorrect(hf,vf,hspotcoff1,vspotcoff1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%horizontal grid line optimizing
[hspotcoff2,hdmean,hdermean]=hgridadjust(hspotcoff1,0,0,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%vertial grid line optimizing
[vspotcoff2,vdmean,vdermean]=vgridadjust(vspotcoff1,0,0,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
et1=cputime-st;
disp(['time after optimizing',num2str(et1)]);
%darw the grid line
[s,hf2]=size(hspotcoff2);
[s,vf2]=size(vspotcoff2);
disp(['draw grid line: ','coordinate y��',num2str(hf2),',coordinate x��',num2str(vf2)]);
%Drawgrid(p,h,w,hspotcoff2,vspotcoff2,hf2,vf2);
%output parameter
para=[et,et1,hf,vf,hf1,vf1,hf2,vf2,hdmean,vdmean,hdermean,vdermean,correct];